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I. Introduction. 

Numerous random number generators have been proposed for the 
System/360. Several of these generators have been incorporated into 
the subroutine library here at the Computer Center. The adequacy of 
some of these generators has rested on the results of some rather weak 
tests for randomness; recent results in the literature have shown many 
of these generators to be very poor performers. This report will 
describe an interim version of a package for random number generation 
which has stood up under intensive statistical testing and is deemed 
to be very satisfactory for the System/760. (The statistical testing 
will be reported elsewhere.) 

The package, when called by a FORTRAN IV program, will deliver 
either a single value or an array (of specified size) of single preci- 
sion, uniforuwly, normally, or exponentially distributed pseudo-random 
deviates or a single value or an array of pseudo-random integers 
uniformly distributed between 1 and go el. The package also has the 
ability (optional) to "shuffle" the pseudo-random nunbers to obtain 
"better" statistical properties. 

Further refinements will be made to this generator; however, it 
is now available for use in an interim version under the name LLRANDOM. 
Future vei«ions will be announced through the W. R. Church Computer 
Center Newsletter. The changes envisioned will be internal and aimed 
at increasing speed and efficiency of coding. The actual numbers pro- 
duced in future versions will remain the same as described here, as will 
the FORTRAN calling sequence. 


Some definitions. By "random number generator" or "pseudo-random 


number generator" is meant an algorithm by which sequences of numbers 


are produced which follow a given probability distribution and possess 
the appearance of randomness. Without attempting to address the still 
unresolved philcsophical question of what a random sequence is, the 
underlined words above are the keys to random number generation on a 
digital computer. The term sequence implies that many random numbers 
must be produced serially from the algorithm. The user may need only a 
very few of these numbers, however we generally require that the algo- 
rithm be able to produce very many numbers. Distribution implies that 
we can associate a probability statement with the occurrence of each 
random number. The distribution is usually taken to be uniform, that 

; is, within a given range the probability of occurrence of a given number 
is the same as for any other number in a similar range. If the algorithm 


produces, say, m distinct numbers then the probability of occurrence 


sa core ta 


for any one of them is 1/m. 

Lastly, we speak of the appearance of randomness. As will be 
shown next, the actual implementation of the algorithm is a recurrence 
relation where each succeeding number is a function of the preceding 
number. True randomness would require independence of successive numbers; 
q however, the algorithm generates a deterministic sequence. Algorithms 

for random number generation do, however, vield sequences which appear 
to be random, hence the term "pseudo-random numbers." It is this 
characteristic which is the subject of statistical testing, that is, 
one asks, "how random does the given sequence appear?" 


The uniform random number generator which forms the basis of 


the :ackage described here is a Lehmer congruential generator. The 


recurrence relation is given by 
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X= A+X | + C (mod m). (1) 


This generator produces integer random numbers between 1 and a. 


These integer values may then be transformed into real-valued numbers 
between 0.06 and 1.0 or into any desired distribution by an appropriate 
: transformation. 
II. The Generator. 

The recurrence relation given in equation (1) is actually called 
f a "Lehmer mixed congruential generator." The term mixed comes from the 
fact that it invc.ves a multiplication by a constant, A, plus an 
addition of a constant, C. The actual implementation used in LLRANDOM 
is called a "multiplicative," or "pure," congruential generator in that 


we take C= 0, giving 


x ZA "Xa (mod m). (2) 


The field of positive integers is, of course, infinite. It is 
a reality of digital computers that only a finite number of positive 
integers are expressible. Specifically, we are limited to the word size : 


of the System/360. This word size is 32 bits with one bit reserved for 


the algebraic sign; hence, an obvious choice for m is 231, The 


product A» xX - is formed by the System/360 in two adjacent registers 


1 4 
yielding a result which may be as large as 263 | We must, however, reduce i F 


this product to a number less than or equal to xe The mod, or modulo, 


operation accomplishes this. The product A °Xo is divided by get 


31 


leaving a quotient which is some multipie of 2°” and a remainder which 


is strictly less than 2s, It is this remainder which is the next 


pseudo-random number Xx in the equation (1). 
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On first examination it would appear that a full 231 numbers 


could be generated by the sequence (1). This is not the case, unless 
A and m are chosen properly. We define a term called the period 
which is the number of unique random numbers computable for a given 
choice of A and m. To illustrate the concept, assume we have a 
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six-bit word with one bit for a sign. We then have m= 2° #= 32. 


Choose A= 9 and work through a sequence starting with Xo = 1. 


Step n X-1 Ae xX -l As X 1 (mod 2°) 
1 1 9 9 
2 9 81 17 
3 17 153 25 
4 25 225 1 
5 1 9 9 


Note that the modulus of this generator is 32, however we have realized 
a period of only 4, that is the sequence of 1, 9, 17, 25,1, 9, . 
repeats after only 4 numbers. U>viously, care must be taken to insure 
that such occurrences do not happen in a random number generator. 
Hopefully, the period will also be independent of the starting value, 
Xo: 
A great deal of work has been done on number theoretic cor idera- 
tions for the choice of m so as to yield a maximum period lengtt 
(see Knuth”), To summarize, generators with modulus m= 2? for ny 
integer, p, can have a maximum period of m/4, or, for the Systeu/360, 


23444 = 229, the period may also depend or the starting value. When 


the modulus m is a prime number, the maximum possible period is m- 1. 
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It so happens that the largest prime less than or equal to 23! is 


31 
2°" - 1, which is most fortuitous. Hence, choosing o = g3t -1 we 


can achieve a maximum period of m-1 =& 231 - 2. These results produce 


only upper bounds on the period length. Recall in the example above, q 

the maximum period possible is 2°/4 = oo 8, but that a period of 

only 4 was observed. This naturally leads to considerations of the ; 

choice of the multiplier, A. ? 
Success in achieving a maximum period lies with the cholce of 

the multiplier. Again, to briefly summarize the pertinent number theory, 

for a modulus 231 the multiplier A must differ by 3 from the nearest 

multiple of 8; the starting value, Xq» must be odd; A _ must be one 

greater than a multiple of 4; and C must be odd. These conditions j 

only assure a maximum period of m/4, not necessarily good statistical ; 

properties. For the random number generator described here (LLRANDOM) i q 

we are choosing C= 0; hence, this length is not achievable if 

m = at Luckily, the conditions on choosing A for the modulus a 

nm = g7t - 1 are more easily met and we can achieve the maximum period. 
Utilizing some of the nice number theoretic properties of the 


ahs 1, to achieve a maximum period, A must be a positive a4 


primitive root of PEL - 1 or a power of such a number. This is generally 4 


number 2 


not easy to find; the value of A used in the generator described here 4 
is -, The number 7 is a positive primitive root of 23} - 1 and raising i 
7 to the fifth power results in the multiplier 16807 which is also a 2 
positive primitive root of ae - 1 (Lewis, Goodman, and Miller‘) and 


satisfies some conditions regarding the statistical performance of the 


generated sequence. These conditions vill not be discussed here. 


The generator 


as 31 
x, = 7 x -1 (mod 2~~-1) (3) 


is the generator reported in Lewis, Goodman, and Miller!. The authors 


cite the results of very extensive tests on this generator, all of which 


show that it is very satisfactory. 
A. Division simulation. A practical consideration for random number 
generators is that they be fast, hopefully without requiring excessive 
Memory to achieve speed. In many applications rather large quartities 
of numbers are needed and the speed of the generator can be crucial. 
Nearly all random number generators are coded as subroutine or 
function subprograms in the assembler or machine language of the computer. 
The algorithm for implementing (3) is rather cimple, involving a multi- 
plication and then a division to effect the modulo operation. On most 
computers the division operation is rather slow as compared to the 
multiplication operation. In the past, the multiplier A was chosen 
so that its binary representation contained many zeroes, thereby speeding 
the multiplication. Unfortunately, this choice was at the expense of 
the period length, since such multipliers rarely met the number theoretic 
conditions for a maximum period. For the LLRANDOM generator (3) described 
here, the division operation has been replaced by a division simulation 
involving two shifts and an add instruction. Should a fixed-point 
(integer) overflow occur, two more additions are required to correct 
the situation. 
The ordinary division on a System/360 Model 67-2 requires 8.49 


micro-seconds. Without overflow, the simulation requires only 3.45 
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micro-seconds. When overflow occurs, the simulation takes an additional 


2.32 micro-seconds for a total of 5.77 micro-seconds. These overflows 


occur quite rarely, on the order of oniy once in 250,000 iterations. 
The division simulation algorithm (again due to Lehmer) is dis- 
cussed by Payne, Rabung, and Bogyo- and works as follows. Define a 


congruence relationship by 
, = e 31 
x’ zA X -l (mod 27"). (4) 


Performing the modulo operation on the product AX 4 would give 


a1 
AX i q2° +4, (5) 


4 where q is some quotient and r is the remainder and is strictly less 


than Neo Adding q to both sides of (4) we get 

oe Eo ges ER Ged 22141) (6) 
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This form gives the desired modulus of 231 - 1, if there is no overflow 


in the addition of x + q. If there Js overflow, to correct it we 


merely add a cunstant of 1 to get 
X = x + q+ 1 (mod 231) =A X-1 (mod a aiy, (7) 


which is again, the desired result. 
This division simulation algorithm is very easily implemented 
on the System/360 and saves considerable execution time over conventional 
division. 
q B. Shuffling. The sequence produced by the generator (3) does appear 


a to consist of independent, uniformly distributed numbers for most purposes. 


We realize that the numbers are not actually independent, due to the 
procedure used to generate *’nem. It has been proposed that a sequence 
of such numbers be further randomized, or "shuffled," to improve upon 
the appearance of randomness (see, for instance, Knuth?) . Serial 
correlation tests are usually employed to detect lack of independence 
in a sequence and at least one generator, RANDU, known to perform badly 
in a three-dimensional serial test was improved by shuffling. These 
tests will be discussed elsewhere. The various shuffling procedures 
which have been put forward have had little empirical validation. 

The package described here has a built-in shuffling mechanism 
and it works as follows. A table of 128 random integers is maintained 
in the package. The starting values in the caple r > resent members of 
the sequence (3) lagged by one million integers starting with an arbi- 
trary seed. When a new integer is generated by the algorithm, its 
right-most seven bits are masked-off to form an index into the table 
(270128). The integer in the table indexed by the right-most seven bits 
is returned to the caller and that table entry is replaced by the 
integer just generated. In essence, we are taking "chunks" of 128 
numbers from the basic sequence and shuffling them before they are used. 

This particular shuffling scheme is dependent on the choice of 
the modulus. For a modulus of ee the right-most bits of a congruential 
random number generator are non-random and their use in this scheme would 
defeat the purpose of shuffling. However, with a modulus of git -1 


and the positive primitive root multiplier A = 7, the right-most 


bits are quite random and the desired results are obtained. 


C. Uniform (0.0,1.0) random numbers. So far, we have discussed how to 


generate uniform random integers over the range 1 to om # 231 -1. In 


most applications, uniform random numbers over the range 0.0 to 1.0 


are desired. In theory, the uniform integers, X are divided by 


ee 


m to produce these numbers, as 


Ue X,/m. (8) 


In actual implementation on the System/360, the integer result is alge- 
braically shifted right seven bits and a normalized floating point expo- 
nent is logically OR'ed on to it. The result is a properly normalized 
floating point random number over the range 0.0 to 1.0, usually referred 
to as a "real" uniform number. 
D. Normall: distributed random deviates. ‘The uniformly distributed 
random numbers described above are not only useful in their own right, 
but form the basis of transformations into random numbers with other 
probability distributions. One of the most important of these distri- 
butions is the Normal distribution. 

There are several methods of approximating a Normal distribution 
with uniform random numbers. One of the oldest and, unfortunately, 
most common is the "sum of k uniforms method." The algorithm is based 
on the fact that the uniform (0.0,1.0) distribution has a mean of 1/2 


and a standard deviation of Y1/12. The algorithm works as follows: 


U, - k/2 
Oe Cities 
Ae/12.0 


X 
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The random deviate X is approximately normally distributed witn mean 
QO and variance 1. The approximation is not as good as other methods 
and it is rather time coneuming in that k uniforms must be generated 
and then summed. It was basically devised to overcome the very time 
cor.suming multiply and divide operations in older computers. 


A more accurate algorithm is known as the Box-Muller method or 


Polar method which is actually a rejection method due to von Neumann. 


The methc1 requires the generation of two uniforms to produce two 
independent Normals. It is based on the distribution of points inside 
the unit circle. The method is more accurate than the "sum of k 
uniforms method" (in fact, theoretically perfect). However, it does 
require two square roots and two natural logarithm operations which are 
generally rather time consuming. 

The algorithm used in the package described here is based on a 
method developed by Marsaglia and is known as the "rectangle-wedge-tail" 
method. This algorithm is by far the fastest algorithm available for 
generating normal]v distributed random numbers, although it requires 
mere memory than the Polar method. 

The second volume of Knuth's “The Art of Computer Programing"” 
gives a complete and detailed description of the algorithm. Briefly, 
the positive half of the Normal density curve is discretized into 37 
rectangles, wedges, and a tail as in Figure 1. All of the rectangles 
are uniformally distributed densities. The wedges are approximated by 
"nearly linear densities." Finally, the tail distribution is computed 


by a modification to the Polar method. The normal density, f(x), is 


then given by the composite function. 
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FIGURE 1 
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f, -f., = "large rectangles" 
= "small rectangles" 
= "wedges" 


f = "tail" 


RECTANGLE-WEDGE-TALL METHOD OF APPROXIMATING THE NORMAL DENSITY 
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f(x) = P,f, Od) + Pof, (x) Ee ctetct P37f47(%), (10) 
wLere 
. 
Pp, = 1, 
i=. } 
the densities f) to £54 are the rectangles; fos to £36 are the 
wedges; and f,. is the tail. The first twelve uniformly distributed 


37 


rectangles are used 88% of the time. This makes for an extremely fast 
algorithm for the majority of deviates. When the taii is sampled, the 
deviate is generated by a modified Polar method and still quite satis- 
factory. 

This generator for Normal deviates, like nearly all others, 
produces deviates with zero mean and unit variance. To change the 
scale and shape to any mean, yy, and the standard deviation, oa, we 


apply the linear transformation 
Z = yp + oX (11) 


where z_ now has the desired shape and scale parameters. 


E. Exponential distribution. Another probability distribution of major 


interest in simulations is the exponential. The cumulative distribu- 
tion function and probability density function for the exponential are 


respectively 
Py) = oer (12) 
E(x) = re ™, (13) 


The expected value of the exponential distriburion is: 


E(X] = 1/a. (14) 
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The problem of generating exponential deviates reduces to one of 
generating "unit" exponentials, i.e. those with A= 1, and then 
multiplying the result by whichever \ is necessary to give the desiced 
distribution. 

One 7c the most common methods of generating numbers from dis- 
tributions zt than the uniform is to use the inverse transformation 


technique (see Gaver and Thompaan’). This can be described graphically, 


as in Figure 2, with a plot of the distribution function 


FIGURE 2. i: 


CUMULATIVE DISTRIBUTION FUNCTION OF THE EXPONENTIAL DISTRIBUTION 


The range of the abscissa, X, is infinite in extent. However, the 
¥ range of the ordinate, F(x), is (0.0,1.0), the range of uniform 
(0,1) random variables. The inverse transformation technique is to 


generate a uniform random number, say U, and use this as the ordinate. 


The exponential deviate, say X, 1s the abscissa point corresponding to 


the intersection of the ordinate and the curve. 


Mathematically, this technique is expressed as 


oeFrGyei- se”, j=l, (15) 
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where u is the uniform random number. This inverse transformation 


is rather easily implemented for exponentially distributed random 


variables via natural logarithms since we get 


pty) = -gn(1-U), 


or by the symmetry of the uniform distribution 
-1 
F(U) = -£n(U)- 


Perhaps the most common impleme atation of exponential deviate 


generators 4s this natural logarithm transformation. It is mathemati- 


cally appealing as well as trivial to program, given the usual FORTRAN 


subroutines. 


The exponential deviate generator in the LLRANDOM package is 


based on Marsaglia's method of dividing the probability density into a 


series of rectangles, wedges, and a tail. Al ‘hough more complicated 


to program and larger in size, this method is approximately 40% faster 


than the logarithmic transformation. 


For a survey of the generation of normal and exponentially 


distributed variables see Ahrens and Dieter . 


III. HOW TO USE THE PACKAGE. 


The random number package described here is intended solely for 


use on the IBM Syscem/360 or System/370 computers. The package consists 


of one Assembler F control section (°SECT) with nine entry points and 


two FORTRAN iV function subprograms. The names of the entry points 


and their functions are summarized in Table 1. 


The subroutine entry point OVFLOW has no calling arguments and 


should be called once and only once at tne beginning of the user's main 


FORTRAN program. The function subprograms RNORTH and REXPTH are called 


by the Assembler routine as needed and should not be called by the 


user. The eight additional entry points are the names of the actual 


routines to generate the random numbers. There are four types of random 


numbers which can be generated: 


(1) uniformly distributed integers on the range 1 to 231 - 1; 


(2) uniformly distributed single precision floating point numbers 
between 0.0 and 1.0; 


(3) single precision floating point normal deviates with mean zero 
and variance 1; and 


(4) single precision floating point exponential deviates with 
mean 1. 


There is a separate entry point for each of the four types if shuffling 


of the sequence is desired. ~ | 


ae ‘ 


For all eight entry points, the FORTRAN calling sequence is the 
same, namely: 

CALL (entry point) (IX, A, N) 
where 

(entry point) refers to the routine desired, 


viz. INT,SINT,RANDOM,SRAND,NORMAL,SNORM,EXPON, or SEXPON; 4 
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is the starting value of the sequence and may contain any integer 
number between 1 and 2147483647. This variable should not be 
altered by the user during the execution of the program, unless 
it is desired to repeat a sequence of random numbers. 


is either a scalar or vector variavle and is the location with 
a specified dimension into which the random number or numbers 
are stored (see next parameter). Note that for entry points 
INT and SINT, this argument should be of INTEGER type. 


is an integer variable or constant designating how many random 
numbers are to be generated during this call. If N is greater 
than 1, A above must be a vector dimensioned at least as large 
as N. If N is equal to 1, then A may be scalar. 


Some sample programs are given below: 


(1) 


To generate 1000 consecutive integer random numbers: 
INTEGER*4 M(1000) 
CALL OVFLOW 


IX = 1234567 


CALL INT (IX, M, 1000) 


END 


To generate 25 shuffled single precision floating point normal 


deviates and scale to mean 10 and standard deviation 5: 
REAL*4 A(25) 

CALL OVFLOW 

JJ = 1936748 


N = 25 


CALL SNORM (JJ, A, N) 
DO 1 I = 1,25 
A(I) = A(I)*5.0 + 10.0 


1 CONTINUE 


END 

(3) To generate one single precision floating point exponential 
deviate with mean 6: 
CALL OVFLOW 


I9 = 98367221 


CALL EXPON (19, E, 1) 


: E = E*6.0 


END 
. A. Implementation. LLRANDOM was designed and coded to run under Operating 
4 System/360 (OS). The Assembler Language control section contains a SPIE 
(Set Program Interrupt Exit) macro instruction which is a part of the OS 


Supervisor Services. This macro enables LLRANDOM to correct for the 


fixed point overflows resulting from the division simulation algorithm. 


The remaincer of the assembly coding is in Basic Assembler Language (BAL), 
i.e. no other macro calls or supervisor calls. To run LLRANDOM under 
another operating system for the System/360, an appropriate substitution 
for the SPIE macro would be necessary. 

As currently programmed, LLRANDOM has the following memory 


requirements: 


SIZE IN BYTES 


MODULE (DECIMAL) 


Assembler CSECT 3571 
FORTRAN function RNORTH isi2 
FORTRAN function REXPTH _1106_ 
Total memory requirement _6189 | 


The System/360 internal timer is rather crude for timing the 
execution of programs. The following times are therefore approximate 
timings for the generation of pseudo-random numbers on a System/360 


Model 67-2. 


ENTRY POINT TIME IN MICROSECONDS 
INT 10.7 
SINT 13.7 
RANDOM 15.6 
SRAND 20.0 
NORMAL 57.5 (Polar method takes 


349 microseconds) 
SNORM 65.8 


EXPON 59.1 (Logarithm method takes 
132 microseconds) 


SEXPON 
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B. Future Enhancements. The normal and exponential deviate routines 


in LLRANDOM are patterned after a package, SUPER-DUPER, available from 
Professor G. Marsaglia at McGill University in Montreal. It is av ilable 
at the Naval Postgraduate School. Marsaglia uses a different multiplier, 
A, and modulus, m, in his congruential generator from that used in 
LLRANDOM and he then exclusive OR's this result with the output of a 
feedback shift register generator. SUPER-DUPER provides onl: one deviate 
per call and does not provide for shuffling the sequence. 

The two FORTRAN function subprograms, RNORTH and REXPTH, are 
taken (with slight modification) directly from SUPER-DUPER. Among the 
changes to be made to LLRANDOM will be to rewrite RNORTH and REXPTH in 
System/360 Assembly Language and incorporate them directly into the 
package. 

We have experienced occasions where la-ge-scale simulations have 
been coded in FORTRAN using double precision variables. The fact that 
LLRANDOM returns single precision numbers causes some inconvenience. 

To alleviate this problen we will provide the capability in LLRANDOM to 
return single precision numbers into double precision variables or arrays. 
Note that the values returned will still be single precision; however, 
they will be stored properly into double precision locations. 

Finally, additional entry points will be added to provide single 
precision floating poin: gamma deviates. Shuffling c.f the gamma deviates 


will also be available. 


Other enhancements are under consideration. 
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